#include "plant2D_SPDE_input_1Dsin.h"

// --------------------------------------------------------------
plant2D_SPDE_input_1Dsin::plant2D_SPDE_input_1Dsin(int M, int L,
						   double idt,double iXYmax){
  int i,m,n;
  // MODEL PARAMETERS -------------------------
  XYmax = iXYmax;

  // NUMERICAL SOLUTION PARAMETERS -----------
  mode = M;
  LatticeSize = L;

  // BASIC BLOCK ------------------------------
  // Input#0 = T (K)
  // Input#1 = W (nm/s)
  // Input#2 = A (nm/s)
  // Input#3 = Freq
  int idim[4] = {1,1,1,1};
  // Output#0 = h (nm/s)
  // Output#1 = r2 (nm^2)
  // Output#2 = m2
  // Output#3 = a (site)
  // Output#4~#7 = values of eigenfunction phi_1 ~ phi_4
  int odim[8];
  odim[0] = 1;
  odim[1] = 1;
  odim[2] = 1;
  odim[3] = 1;
  odim[4] = (mode+1)*(mode+1); // The mode with sin*sin
  odim[5] = (mode+1)*(mode+1);
  odim[6] = (mode+1)*(mode+1);
  odim[7] = (mode+1)*(mode+1);
  block::initialize(4,idim,8,odim,idt);

  // STATES -----------------------------------
  z1 = block::outputs[4];
  z2 = block::outputs[5];
  z3 = block::outputs[6];
  z4 = block::outputs[7];
  z1_next = new double[(mode+1)*(mode+1)];
  z2_next = new double[(mode+1)*(mode+1)];
  z3_next = new double[(mode+1)*(mode+1)];
  z4_next = new double[(mode+1)*(mode+1)];
  for(i=0;i<(mode+1)*(mode+1);i++){
    z1_next[i] = 0.0;
    z2_next[i] = 0.0;
    z3_next[i] = 0.0;
    z4_next[i] = 0.0;
  }

  // MODEL PARAMETERS -----------------------------------
  lambda = new double*[mode+1];
  for(m=0;m<mode+1;m++){
    lambda[m] = new double[mode+1];
  }

  for(i=0;i<4;i++){
    Km2[i] = new double[(mode+1)*(mode+1)];
  }

  double dphi;
  int ix,iy,ix_next;
  int index;
  double temp;
  for(m=0;m<=mode;m++){
    temp = 4*sin(pi*m/LatticeSize)*sin(pi*m/LatticeSize)/(XYmax*XYmax);
    for(n=0;n<=mode;n++){
      index = index2D(m,n);
      //Version 1
      //for(i=0;i<4;i++){
      //  Km2[i][index] = 0.0;
      //}
      //for(ix=0;ix<LatticeSize;ix++){
      //  for(iy=0;iy<LatticeSize;iy++){
      //    if(ix == LatticeSize-1){ix_next = 0;}
      //    else{ix_next = ix+1;}
      //    dphi = phi1mn(m,n,ix_next,iy)-phi1mn(m,n,ix,iy);
      //    Km2[0][index] += dphi*dphi;
      //    dphi = phi2mn(m,n,ix_next,iy)-phi2mn(m,n,ix,iy);
      //    Km2[1][index] += dphi*dphi;
      //    dphi = phi3mn(m,n,ix_next,iy)-phi3mn(m,n,ix,iy);
      //    Km2[2][index] += dphi*dphi;
      //    dphi = phi4mn(m,n,ix_next,iy)-phi4mn(m,n,ix,iy);
      //    Km2[3][index] += dphi*dphi;
      //  }
      //}
      //Km2[0][index] = Km2[0][index]/(XYmax*XYmax); //FIXME
      //Km2[1][index] = Km2[1][index]/(XYmax*XYmax);
      //Km2[2][index] = Km2[2][index]/(XYmax*XYmax);
      //Km2[3][index] = Km2[3][index]/(XYmax*XYmax);

      // Version 2
      //Km2[0][index] = 4*sin(pi*m/LatticeSize)*sin(pi*m/LatticeSize)/(XYmax*XYmax);
      //Km2[1][index] = Km2[0][index];
      //Km2[2][index] = Km2[0][index];
      //Km2[3][index] = Km2[0][index];
      
      // Version 3

      for(i=0;i<4;i++){
	switch(i){
	case 0:
	  if(m==0||2*m==LatticeSize||n==0||2*n==LatticeSize){
	    Km2[i][index] = 0.0;
	  }else{
	    Km2[i][index] = temp;
	  }
	  break;
	case 1:
	  if(m==0){
	    Km2[i][index] = 0.0;
	  }else{
	    Km2[i][index] = temp;
	  }
	  break;
	case 2:
	  if(m==0||2*m==LatticeSize){
	    Km2[i][index] = 0.0;
	  }else{
	    Km2[i][index] = temp;
	  }
	  break;
	case 3:
	  if(n==0||2*n==LatticeSize||m==0){
	    Km2[i][index] = 0.0;
	  }else{
	    Km2[i][index] = temp;
	  }
	}
      }
    }
  }

  // Simulatin Parameters -------------------------
  // For Analytical Expression
  sqrt_dt = sqrt(block::dt);
  // For explicit/implicit method
  dt_con = 0.0001;
  sqrt_dt_con = sqrt(dt_con);
  dt_conservative_factor = 0.1;
  
  gsl_rng_env_setup();
  noise_type = gsl_rng_default;
  noise = gsl_rng_alloc(noise_type);
  gsl_rng_set(noise,time_seed());
  
  // INITIALIZATION -------------------------------
  reset();
}
// --------------------------------------------------------------
plant2D_SPDE_input_1Dsin::~plant2D_SPDE_input_1Dsin(){
  int m;
  for(int i=0;i<4;i++){
    delete []Km2[i];
  }

  for(m=0;m<mode+1;m++){
    delete []lambda[m];
  }

  delete []z1_next;
  delete []z2_next;
  delete []z3_next;
  delete []z4_next;

  delete []lambda;
  gsl_rng_free(noise);
}
// --------------------------------------------------------------
void plant2D_SPDE_input_1Dsin::reset(){
  block::time = 0.0;
  int m,n;
  // Reset states
  for(m=0;m<mode+1;m++){
    for(n=0;n<mode+1;n++){
      z1[index2D(m,n)] = 0.0;
      z2[index2D(m,n)] = 0.0;
      z3[index2D(m,n)] = 0.0;
      z4[index2D(m,n)] = 0.0;
    }
  }
  // Reset outputs
  outputs[0][0] = getH();
  outputs[1][0] = getR2();
  outputs[2][0] = getM2();
  outputs[3][0] = getCorLength();
}
// --------------------------------------------------------------
void plant2D_SPDE_input_1Dsin::update(double sysTime){
  if(sysTime<time){
    cerr << "Error@plant2D_SPDE_input_1Dsin::update(), sysTime = " << sysTime << ", time = " << time << "\n";
    exit(1);
  }

  // Update states  ====================================
  int i,m,n;
  int index;
  double c2;
  double sigma2;
  double iT,iWmean;
 
  while(time+dt<=sysTime || time == 0){
    if(sysTime!=0){
      // x_{k} <== x_{k+1}
      for(i=0;i<(mode+1)*(mode+1);i++){
	z1[i] = z1_next[i];
	z2[i] = z2_next[i];
	z3[i] = z3_next[i];
	z4[i] = z4_next[i];
      } 
    }

    // Calculate x_{k+1} ---------------------------
    iT = inputs[0][0];
    iWmean = inputs[1][0];
    A = inputs[2][0];
    Freq = inputs[3][0];

    if(iT>1500 || iT<273){
      cerr << "Invalid value of T: " << iT << "\n";
      exit(1);
    }

    if(iWmean>10.0 || iWmean<0.0){
      cerr << "Invalid value of W: " << iWmean << "\n";
      exit(1);
    }

    if(A>iWmean && A<0.0){
      cerr << "Invalid value of A: " << A << "\n";
      exit(1);
    }

    if(Freq<0){
      cerr << "Invalid value of Freq: " << Freq << "\n";
      exit(1);
    }
    // UPDATE MODEL PARAMETERS ---------------------
    if(iT!=T || iWmean!=W){
      T = iT;
      W = iWmean;
      c2 = eq_c2(T,W);
      // sigma2 = sqrt(eq_sigma2(T,W)); // A BIG BUG
      sigma2 = eq_sigma2(T,W);
      for(m=0;m<mode+1;m++){
	for(n=0;n<mode+1;n++){
	  lambda[m][n] = -4*c2*(m*m+n*n)*pi*pi/(XYmax*XYmax);
	  //	  dt_con[m][n] = -2/lambda[m][n]*dt_conservative_factor;
	}
      }  
    }
    else{
      c2 = eq_c2(T,W);
      sigma2 = eq_sigma2(T,W);
    }

    // UPDATE STATES -------------------------------
    for(m=0;m<=mode;m++){
      for(n=0;n<=mode;n++){
	index = index2D(m,n);

	// Implicit method =====================================
	/*	if(m>0&&n>0){
	  // z1[index] = z1[index]+lambda[m][n]*z1[index]*dt_con+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con;
	  z1[index] = (z1[index]+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con)/(1-lambda[m][n]*dt_con);
	}

	if(m==0&&n==0){
	  //	  z2[index] = z2[index]+c_200*dt_con+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con;
	  double w200 = iWmean*XYmax*(1-(cos(2*iFreq*pi)-1)/(2*iFreq*pi));
	  z2[index] = (z2[index]+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con+w200)/(1-lambda[m][n]*dt_con);
	}
	else if(n!=0 || iFreq==m){
	  //	  z2[index] = z2[index]+lambda[m][n]*z2[index]*dt_con+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con;
	  z2[index] = (z2[index]+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con)/(1-lambda[m][n]*dt_con);
	}
	else{// n = 0 && iFreq!=m
	  double w2m0 = SQRT_TWO*iWmean*XYmax*iFreq*cos(2*iFreq*pi)/(2*pi*(m*m-iFreq*iFreq));
	  z2[index] = (z2[index]+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con+w2m0)/(1-lambda[m][n]*dt_con);
	}
     
	if(m>0){
	  if(n!=0){
	    //	    z3[index] = z3[index]+lambda[m][n]*z3[index]*dt_con+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con;
	    z3[index] = (z3[index]+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con)/(1-lambda[m][n]*dt_con);
	  }
	  else{
	    if(m==iFreq){
	      double w3k0 = SQRT_TWO*iWmean*XYmax/2;
	      z3[index] = (z3[index]+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con+w3k0)/(1-lambda[m][n]*dt_con);
	    }
	    else{//m!=iFreq
	      double w3m0 = SQRT_TWO*iWmean*XYmax*m*sin(2*iFreq*pi)/(2*pi*(iFreq*iFreq-m*m));
	      z3[index] = (z3[index]+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con+w3m0)/(1-lambda[m][n]*dt_con);
	    }
	  }
	 
	  If(n>0){
	    //	  z4[index] = z4[index]+lambda[m][n]*z4[index2D(m,n)]*dt_con+gsl_ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con;
	    z4[index] = (z4[index]+gsl-ran_gaussian_ziggurat(noise,1.0)*sigma2*sqrt_dt_con)/(1-lambda[m][n]*dt_con);
	  }
	  }*/
	// Analytical expression ==============================
	double temp1 = exp(lambda[m][n]*dt);
	double temp2 = sqrt(sigma2*(exp(2.0*lambda[m][n]*dt)-1.0)/(2.0*lambda[m][n]));
	double temp3 = (temp1-1.0)/lambda[m][n];
	
	if(m>0&&n>0){
	  z1_next[index] = z1[index]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
	}
	
	if(m==0&&n==0){
	  // double w200 = W*XYmax+A*XYmax*(1-cos(2*Freq*pi))/(2*Freq*pi);
	  double w200 = W*XYmax;
	  z2_next[index] = z2[index]+w200*dt+gsl_ran_gaussian_ziggurat(noise,1.0)*sqrt(sigma2)*sqrt_dt;
	}
	else{
	  if(n==0 && Freq!=m){
	    // double w2m0 = SQRT_TWO*A*XYmax*Freq*(cos(2*Freq*pi)-1)/(2*pi*(m*m-Freq*Freq));
	    double w2m0 = 0.0;
	    z2_next[index] = z2[index]*temp1+
	      temp2*gsl_ran_gaussian_ziggurat(noise,1.0)+temp3*w2m0;
	  }
	  else{
	    z2_next[index] = z2[index]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
	  }
	}

	if(m>0){
	  if(n!=0){
	    z3_next[index] = z3[index]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
	  }
	  else{
	    if(m==Freq){
	      double w3k0 = SQRT_TWO*A*XYmax/2;
	      z3_next[index] = z3[index]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0)+temp3*w3k0;
	    }
	    else{
	      // double w3m0 = SQRT_TWO*A*XYmax*m*sin(2*Freq*pi)/(2*pi*(Freq*Freq-m*m));
	      double w3m0 = 0.0;
	      z3_next[index] = z3[index]*temp1+
		temp2*gsl_ran_gaussian_ziggurat(noise,1.0)+temp3*w3m0;
	    }
	  }
	}
	else{
	  z3_next[index] = 0.0;
	}
	
	if(n>0){
	  z4_next[index] = z4[index]*temp1+temp2*gsl_ran_gaussian_ziggurat(noise,1.0);
	}
	else{
	  z4_next[index] = 0.0;
	}
      } // end of for(n= ...)
    }// end of for(m = ...)
    // ==========================================
    if(time==0)
      break;
    else{
      time=time+dt;
    }
    // ==========================================
  }// end of while

  // OUTPUTS -----------------------------------
  outputs[0][0] = getH();
  outputs[1][0] = getR2();
  outputs[2][0] = getM2();
  outputs[3][0] = getCorLength();
  return;
}
// --------------------------------------------------------------
void plant2D_SPDE_input_1Dsin::output(string FileName){
  ofstream fout;
  fout.open(FileName.c_str());
  fout << "Mode = " << mode;
  fout << "Time = " << block::time;
  int i;
  
  int dim = (mode+1)*(mode+1);
  
  for(i=0;i<dim;i++){
    fout << z1[i] << " ";
  }
  fout << "\n";
  
  for(i=0;i<dim;i++){
    fout << z2[i] << " ";
  }
  fout << "\n";

  for(i=0;i<dim;i++){
    fout << z3[i] << " ";
  }
  fout << "\n";
  
  for(i=0;i<dim;i++){
    fout << z4[i] << " ";
  }
}
// --------------------------------------------------------------
double plant2D_SPDE_input_1Dsin::getM2(){
  int m,n,index;
  double m2 = 0.0;
  for(m=0;m<=mode;m++){
    for(n=0;n<=mode;n++){
      index = index2D(m,n);
      m2 +=  Km2[0][index]*z1[index]*z1[index]
	+Km2[1][index]*z2[index]*z2[index]
	+Km2[2][index]*z3[index]*z3[index]
	+Km2[3][index]*z4[index]*z4[index];
    }
  }
  return m2;
}
// --------------------------------------------------------------
double plant2D_SPDE_input_1Dsin::getR2(){
  int m,n,index;
  double r2 = 0.0;
  for(m=0;m<=mode;m++){
    for(n=0;n<=mode;n++){
      if(m!=0||n!=0){
	index = index2D(m,n);
	r2 +=  z1[index]*z1[index]
	  +z2[index]*z2[index]
	  +z3[index]*z3[index]
	  +z4[index]*z4[index];
      }
    }
  }
  r2 /= (XYmax*XYmax);
  return r2;
}
// --------------------------------------------------------------
double plant2D_SPDE_input_1Dsin::getH(){
  double hmean;
  hmean = z2[index2D(0,0)]/XYmax;
  return hmean;
}
// --------------------------------------------------------------
double plant2D_SPDE_input_1Dsin::getCorLength(){
  double r2 = outputs[1][0];
  double m2 = outputs[2][0];
  return sqrt(2*r2/m2);
}
// --------------------------------------------------------------
void plant2D_SPDE_input_1Dsin::getSnapshot(string FileName){
  ofstream fsnapshot(FileName.c_str());
  int ix,iy,im,in;
  double h;
  int index;
  for(ix=0;ix<LatticeSize;ix++){
    for(iy=0;iy<LatticeSize;iy++){
      h = 0;
      for(im=0;im<=mode;im++){
	for(in=0;in<=mode;in++){
	  index = index2D(im,in);
	  h += 
	    phi1mn(im,in,ix,iy)*z1[index]+
	    phi2mn(im,in,ix,iy)*z2[index]+
	    phi3mn(im,in,ix,iy)*z3[index]+
	    phi4mn(im,in,ix,iy)*z4[index];
	}
      }
      fsnapshot << setw(10) << h;
    }
    if(ix!=LatticeSize-1){
      fsnapshot << "\n";
    }
  }
  fsnapshot.close();
}
// --------------------------------------------------------------
double plant2D_SPDE_input_1Dsin::eq_c2(double T,double W){
  // INPUTS:
  // T: k
  // W: nm/s
  //
  // OUTPUTS:
  // C2: nm^2/s
  //
  //  double pc1 = 7.456e-8;
  //  double pc2 = -4.039e-7;
  //  double pc3 = 8.092e-7;
  //  double pc4 = -7.376e-7;
  //  double pc5 = 7.579e-7;
  //  return pc1*W*W*W*W+pc2*W*W*W+pc3*W*W+pc4*W+pc5;
  //  return k_c*exp(k_B*T/E_c)*pow(W,-a_c)/(LatticeSize*LatticeSize);
  double pc1 = -10.53;
  double pc2 = 12.23;
  double pc3 = -4.953;
  double pc4 = 0.8804;
  return pc1*W*W*W+pc2*W*W+pc3*W+pc4;
}

double plant2D_SPDE_input_1Dsin::eq_sigma2(double T,double W){
  //  double ps1 = 1.224e-4;
  //  double ps2 = 2.871e-6;
  //  return ps1*W+ps2;

  //  double index = (a_t+k_t*W)*T-a_v-k_v*W;
  //  return pi*pi/(LatticeSize*LatticeSize)*(1+exp(index));

  // 
  double ps1 = -15.98;
  double ps2 = 16.39;
  double ps3 = -5.634;
  double ps4 = 0.8773;
  double ps5 = 0.01836;

  return ps1*W*W*W*W+ps2*W*W*W+ps3*W*W+ps4*W+ps5;
}
// -------------------------------------------------------------
